function [F, G, obj] = FNMTF(X, nFeaClusters, nSmpClusters)
% Fast Nonnegative Matrix Tri-Factorization 
% Input
%         X: nDim * nSmp
%         nFeaClusters
%         nSmpClusters
% Output
%         G: nSmp * nSmpClusters
%         F: nDim * nFeaClusters
% min ||X - F S G'||^2
% s.t. F \in {0,1}^{d,m}; G \in {0,1}^{n, c}; 
% 
% Note
% (1). F, G are cluster indicator matrices, Not relaxed other matrices
% (2). S is not constrained to be Non-negative, which means this algorithms
% can be applied when X has negative elements
% (3). The optimization of F, G are decoupled, which makes it is easy to
% optimization. 
% 
% [1]. Fast Nonnegative Matrix Tri-Factorization for Large-Scale Data
% Co-Clustering, Hua Wang, IJCAI, 2011
% 

[nDim, nSmp] = size(X);

% **********************************************************
% Initialize G and F with arbitrary class indicator matrices
% **********************************************************

r = randi(nFeaClusters, [1, nDim]);
F = zeros(nDim, nFeaClusters);
F(sub2ind([nDim, nFeaClusters], 1:nDim, r')) = 1;

r = randi(nSmpClusters, [1, nSmp]);
G = zeros(nSmp, nSmpClusters);
G(sub2ind([nSmp, nSmpClusters], 1:nSmp, r')) = 1;

nIter = 200;
obj = zero(nIter,1);
epsilon = 1e-5;
for iter = 1:nIter
    
    % calculate S by Eq (6) ;
    S = F' * X * G;
    S = diag(1./sum(F)) * S * diag(1./sum(G));
    
    % calculate G by Eq. (8) ;
    smpCenter = F * S;
    smpIdx = knnsearch(smpCenter', X', 'K', 1);
    G(sub2ind([nSmp, nSmpClusters], 1:nSmp, smpIdx')) = 1;
    
    % calculate F by Eq. (8) ;
    feaCenter = S * G';
    feaIdx = knnsearch(feaCenter, X, 'K', 1);
    F(sub2ind([nDim, nFeaClusters], 1:nDim, feaIdx')) = 1;
    
    % check convergence
    if iter > 1;
        if abs(obj(end) - obj(end-1))/abs(obj(end)) < epsilon;
            obj(iter+1:end) = [];
            break;
        end
    end
end

if iter == nIter
    warning('Not Converged! More iteration needed');
end

return;
    
